using Pkg
Pkg.activate(".")

using DataFrames, CSV, StatsPlots, Downloads, DataFramesMeta, Dates, Statistics
using Turing, Memoization, ReverseDiff, Serialization
using LinearAlgebra, Distributions

Firearms and Violent Crime, recent trends

In late June 2022 the Supreme Court of the United States (SCOTUS) issued a ruling in New York State Rifle and Pistol Association v Bruen which ended the possibility for any state to issue concealed weapons (CCW) permits on the basis of judgement by government employees (known as "may issue" laws). For example in New York an applicant had to convince a police department that they had a "special need" for protection (such as for example a credible death threat) and in California the Sherriff of each county had a similar discretionary ability to decide whether there was a "sufficient cause" for issuing a permit. Simple self defense was not enough.

In the Bruen ruling, either permitless systems as now exist in 25 states, or CCW permit systems in which "objective" standards apply were the only allowable options (known as "shall issue" laws). Objective standards are for example "has never been convicted of a felony" or "has passed a required safety class" and other similar requirements which are factually verifiable by anyone.

After this ruling, as with most rulings liberalizing gun laws, political opponents of the ruling predicted "wild west" style shootouts in the streets. But in fact at the time of the ruling 43 states had either permitless carry (25 states) or a "shall issue" CCW permit system (18 states).

The usconcealedcarry.com website shows the following "shall issue" map. The only misleading data point is that Vermont is shown as red, because they have never regulated concealed carry in any way, it was explicitly constitutionally guaranteed from the beginning and so there is no CCW permit to issue at all.

Shall issue states

So what do we know about firearms carry and violence?

Let's begin by loading several datasets. We have data on background checks each month, data matching FIPS codes for states, and CDC WONDER data on firearms homicides, as well as all homicides. We separate homicides from firearms suicides as they have a different cause and require different prevention techniques than violent crime. In general firearms suicides are about 2x as big as all homicides.

Loading and munging the data into useful form:

if !Base.Filesystem.ispath("data/nics-firearm-background-checks.csv")
    Downloads.download("https://github.com/BuzzFeedNews/nics-firearm-background-checks/raw/master/data/nics-firearm-background-checks.csv","./data/nics-firearm-background-checks.csv")
end

if ~Base.Filesystem.ispath("data/fipscodes.csv")
    Downloads.download("https://www2.census.gov/geo/docs/reference/state.txt","./data/fipscodes.csv")
end
fipscodes = @chain CSV.read("data/fipscodes.csv",DataFrame) begin  
    @subset(:STATE .< 60)
end


# manually downloaded CDC wonder gun homicide data using online form 
gunhomicidesnew = CSV.read("./data/wonder-gun-homicide-byyear.csv",DataFrame)
gunhomicidesold = CSV.read("./data/cdc-wonder-firearm-homicide-1979-1998.csv",DataFrame)
gunhomicides = [gunhomicidesnew ; gunhomicidesold]
gunhomicides = @chain rename(gunhomicides,Dict("Crude Rate" => "crudedeathrate", "State Code" => "StateCode", "Year" => "year", "Year Code" => "YearCode")) begin
    @subset(.! ismissing.(:Deaths))
    @transform(:crudedeathrate = map(x -> isnothing(x) ? missing : x, tryparse.(Float64,String.(:crudedeathrate))),
        :state = String31.(:State))
    @subset(.! ismissing.(:crudedeathrate))
    @orderby(:StateCode,:year)
end


gunchecks = CSV.read("./data/nics-firearm-background-checks.csv",DataFrame)
startpops = @chain gunhomicides, @subset(:year .== 1999), @select(:State, :StartingPop = :Population,:StartDeathRate = :crudedeathrate)



gunsalesest = @chain gunchecks,
    @orderby(:state,:month),
    @select(:state,:month,:salesest = 1.1 .* (:handgun .+ :long_gun) .+ 2.0 .* :multiple),
    groupby(:state),
    @transform(:cumsales = cumsum(:salesest),:year=year.(:month)),
    @subset(month.(:month) .== 6)


joined = leftjoin(gunhomicides,gunsalesest, on=[:State => :state,:year => :year]; matchmissing=:notequal,makeunique=true)

joined = @chain leftjoin(joined,startpops; on=[:state => :State],makeunique=true) begin
    @subset(.! ismissing.(:Deaths) .&& .! ismissing.(:StartingPop))
    @subset(:StateCode .< 60)
    @transform(:cumgunrate = (:cumsales .+ 1.0 .* :StartingPop) ./ :Population )
end

joined = @orderby(leftjoin(joined,fipscodes,on = :state => :STATE_NAME),:StateCode,:year)

Let's first take a look at total guns vs gun homicides

Although we don't have a proper starting value for the beginning of the period, we know that for quite a while now the number of guns per capita in the US is on the order of 1. We therefore put a starting value of 1 for every state even though in fact this should vary somewhat likely between 0.7 and 1.5 or so, and then cumulatively total up the approximate number of firearms purchased since then and divide by the running population estimate, to calculate a firearms per capita time series.

We then plot firearms homicides vs firearms per capita as a parametric curve with time not visible. If purchasing guns causes gun homicides, then as guns per capita increases through time so will firearms homicides. However there are also other models that could cause such trends.

Nevertheless, if gun purchases do cause gun homicides, we should see these positive trending lines.

display(plot(joined.cumgunrate,joined.crudedeathrate; group=joined.state, legend=false, xlab="Gun Rate (guns/capita)",
    ylab="Firearm Homicide Rate (deaths/100kcapita/yr)", title="Gun Homicide vs Gun Prevalence", alpha=.5))

display(plot(joined.cumgunrate,joined.crudedeathrate; group=joined.state, legend=false, xlab="Gun Rate (guns/capita)",
    ylab="Firearm Homicide Rate (deaths/100kcapita/yr)", title="Gun Homicide vs Gun Prevalence", alpha=.5, ylim=(0,12)))

Firearms data as a time series

In addition we plot how Homicides have changed through time, and how firearms per capita has changed through time 

display(plot(joined.year,joined.crudedeathrate; group=joined.state, legend=false, xlab="Year",
    ylab="Firearm Homicide Rate (deaths/100kcapita/yr)", title="Gun Homicide vs Time", alpha=.5))

display(plot(joined.year, joined.cumgunrate; group=joined.state, legend=false, xlab="Year",
    ylab="Firearms Per Capita Estimate (guns/person)", title="Gun Ownership vs Time", alpha=.5, ylim=(0,4),xlim=(1999,2020)))

A zoom in to regional trends

Constitutional Carry

The most recent trend in firearms carry in the US is the conversion of "shall issue" states to a permitless carry system. In this system everyone over a certain age (18 or 21 depending on state) who is legally able to possess a pistol may carry that pistol concealed in public subject to relatively few restrictions (such as not allowed in secure areas with metal detectors etc). For the moment we ignore the issue of "open carry" which has different issues such as intimidation purposes.

Of all the systems for regulating firearms in public, this is the least restrictive. There are various models one might have of how this would affect gun homicides. Some might argue armed citizens deter crime, others might argue it makes it easier to commit crime, others might argue the effect if any would be minimal as criminals were already carrying guns and only relatively "model citizens" would start carrying, and they are generally not involved in crime in terms of committing them or as victims. It's worthwhile to look at the data to see what hypotheses are comopatible with the observed trends.

We plot firearms homicide rates through time for each state. We place a vertical red line at the point in time where the constitutional carry law was signed if any.

concarry = leftjoin(fipscodes,CSV.read("./data/ConstCarryDates.csv",DataFrame),on = :STUSAB => :StateCode)
concarry.date = map(x -> ismissing(x) ? missing : Date(div(x,10000),div(mod(x,10000),100),1),concarry.LawDate)
concarry.year = map(x -> ismissing(x) ? missing : year(x),concarry.date)
concarry = @subset(concarry,:STATE .<= 59) #ignore minor territories

let p = []
    for (st,stusab,statename,statens,lawdate,stdate,year) in eachrow(concarry)
        d = year
        sub = @subset(joined,:state .== statename)
        if(nrow(sub) > 0)
            pp = plot(sub.year,sub.crudedeathrate,ylim=(0,15),size=(500,500),
                title="$statename\nhomicides",xlab="Year",ylab="Gun Homicide/100k/yr",legend=false)
            if(!ismissing(d))
                pp = plot!([(d,0),(d,15)])
            end
        push!(p,pp)
        end
    end
    display(plot(p...; size = (2000,2000)))
end

States by Geographic Group

Looking at the raw data we can see that there are a number of states that have an "accelerating" trend of increased firearms homicide violence in the period from about mid 2014 to 2015 onward. Some of these have converted to Constitutional Carry near the beginning of that period or during that period, with the trend beginning before the carry provision for several of the states. Also many states have passed their laws in the last few years and have no data post-law.

Let's look at these states grouped into geographic regions, because trends in regions may be more obvious. We will start with the following groupings:

stategroups = [["WA","OR","CA","NV","AZ"],["ID","MT","WY","UT","CO","NM","TX"],["ND","SD","NE","KS","OK","IA","MO"],
    ["MN","WI","IL","IN","MI","OH"],["AR","LA","MS","AL","GA","FL","SC"],["TN","NC","KY","WV","VA","MD","DE"],
    ["PA","NJ","NY","CT","RI","MA","VT","NH","ME"],["AK","HI","DC"]]

stategroupdf = DataFrame(STUSAB = reduce(vcat,stategroups),stgroup = reduce(vcat,[[k for j in 1:length(stategroups[k])] for k in 1:length(stategroups)]))
stategroupdf = leftjoin(fipscodes,stategroupdf; on = :STUSAB)

@assert(length(unique(reduce(vcat,stategroups))) == 51) # 50 states plus DC

let p1 = [], p2=[]
    for statelist in stategroups
        sub = joined[in.(joined.STUSAB, Ref(statelist)),:]
        ccdates = @chain @subset(concarry, in.(concarry.STUSAB,Ref(statelist))) begin
            @subset(.! ismissing.(:year))
        end
        
        ccdatesenact = leftjoin(ccdates,sub, on = [:STUSAB => :STUSAB, :year => :year], makeunique = true)
        ccdatesenact = @subset(ccdatesenact, .! ismissing.(:crudedeathrate))
        #display(ccdatesenact)
        #display(sub)
        p = plot(sub.year,sub.crudedeathrate,group=sub.STUSAB,ylim=(0,20),legend=:topleft,linewidth = 3, thickness_scaling=1,
            ylab="Homicides/100k population/yr",xlab="Year")
        if nrow(ccdatesenact) > 0
            p = scatter!(ccdatesenact.year,ccdatesenact.crudedeathrate,markersize=6,group = ccdatesenact.STUSAB)
        end
        push!(p1,p)
        p = plot(sub.cumgunrate,sub.crudedeathrate,markersize=6,linewidth=3,thickness_scaling=1, xlab="Gun Ownership (guns/person)", ylab="Homicide rate (per 100k)", group=sub.STUSAB)
        push!(p2,p)
    end
    display(plot(p1...; size=(1000,1000)))
    display(plot(p2...; size=(1000,1000)))
end

States by Relative Change

Each state has its own overall level of firearm violence, but often trends in the same direction relative to its region even if its overall level is higher or lower. This suggests a model in which we measure each state using a dimensionless number relative to some baseline level. For example the average rate in the years 1999,2000,2001

baselines = @chain joined begin
    @subset(in.(:year,Ref((1999,2000,2001))))
    @by(:STUSAB, :baseline = mean(:crudedeathrate))
end

bljoined = @chain leftjoin(stategroupdf,leftjoin(baselines,joined,on = :STUSAB, makeunique=true),on=:STUSAB,makeunique=true) begin
    @subset(.! ismissing.(:crudedeathrate))
    @orderby(:stgroup,:STUSAB,:year)
end

function getccd(gr)
    ccd = innerjoin(gr,@subset(concarry,.! ismissing.(:year)), on = [:STUSAB => :STUSAB, :year => :year], makeunique = true)
    ccd = @subset(ccd, .! ismissing.(:crudedeathrate))
    ccd
end

function makerelplot(gr)
    ccd = getccd(gr)
#    display(ccd)
    p = plot(gr.year,gr.crudedeathrate ./ gr.baseline,group=gr.STUSAB,legend=:topleft,ylim=(0.0,2.5),xlab="Year",ylab="Relative Rate") 
    if nrow(ccd) == 1
        scatter!(ccd.year,ccd.crudedeathrate ./ ccd.baseline)
    elseif nrow(ccd) > 1 
        scatter!(ccd.year,ccd.crudedeathrate ./ ccd.baseline,group = ccd.STUSAB)
    end
    p
end

makerelplot (generic function with 1 method)

What has been the relative change?

Each state plotted in a group with other geographically similar states. Firearm homicide rate relative to average recorded rate 1999-2001.

plots = [makerelplot(gr)
    for gr in groupby(bljoined,:stgroup)]
    
plot(plots... ; size=(1000,1000),linewidth=2,xlab="Year",ylab="Relative Rate")

Modeling the process with Bayesian Models

We can try to estimate the effects of these laws by building a model of the overall process. Because we can't resurrect people from the dead, the homicide rate can never be negative. It makes sense then to model the homicide rate on a logarithmic scale.

We estimate the behavior of the states as an overall level, plus a shape that is common to the region. Individual states are then allowed a small perturbation to the regional shape.

This small perturbation method was the only way that seemed to make sense to estimate the counterfactual. Clearly none of the states will perform exactly like the average. So if we use the regional average, we will bias towards a lot of noise. Obviously, there is no effect if no law was passed, so for those states with no law passed, the counterfactual is the same as the actual, and we then should allow a size of perturbation which allows the function to fit the actuals for those states. By modifying the parameter statecoefpert we can make it as small as possible while still causing all the non-law states to fit reasonably well. At this point whatever lack of fit there is post-law in those states, could be at least partially attributed to the law. The choice of value of log(1.25)/2.0 ~ 0.112 seemed to fit the bill, leaving only Maryland the only non-law not quite well fit by the model. Any smaller value left both MD and NY fitting poorly.

The effect of constitutional carry law would then be estimated as the log firearm homicide rate which exceeds the state level counterfactual estimate in the years after the passage of the law, with a transient onset curve included in the model.

We impose an informal smoothness requirement on the counterfactual estimates by using a compact radial basis function expansion with one center every 4 years, and a maximum radius of 7 years so that there is overlap between adjacent centers. 

function bumpfun(x,c,scale)
    stdx = (x-c)/scale
    if stdx < -1.0 || stdx > 1.0
        0.0
    else
        exp(1.0-1.0/(1.0-stdx^2)) # goes to zero at -1 and 1, and 1 at x=0
    end
end

function timeser(x,coefs,centers,scale)
    f = 0.0
    for (a,c) in Iterators.zip(coefs,centers)
        f += a*bumpfun(x,c,scale)
    end
    return f
end

function laweffect(yr,rate,start)
    if yr >= start
        1.0-exp(-rate*(yr-start))
    else   
        0.0
    end
end


statecoefpert = log(1.25)/2.0

include("model.jl") ## load the model, this is a separate file so we can save the output unless it changes

modeljoined = @chain leftjoin(joined,stategroupdf; on = :STUSAB, makeunique=true) begin
    @subset(.! ismissing.(:crudedeathrate))
end

lawdates = @chain leftjoin(DataFrame(STATE=1:55),concarry; on=:STATE) begin
    @subset(:STATE .< 60)
    @rtransform(:year = ismissing(:year) ? 3000 : :year)
    @orderby(:STATE)
end

centers = [i for i in (minimum(modeljoined.year)-4):4:(maximum(modeljoined.year)+4)]
width = 7.0


modl = guns(modeljoined.StateCode,maximum(modeljoined.StateCode),modeljoined.stgroup,maximum(modeljoined.stgroup),
        centers, width, modeljoined.year,log.(modeljoined.crudedeathrate),lawdates.year,statecoefpert)

setadbackend(:reversediff)
Turing.setrdcache(true)

savedfile = "./saved/samples.dat"
global s = []
if stat(savedfile).mtime > stat("./model.jl").mtime
    global s = deserialize(savedfile)
else
    global s = sample(modl,NUTS(500,0.8),MCMCThreads(),500,3)
    serialize(savedfile,s)
end

Having sampled the model, we can then compute summary graphs of the results. We will show the posterior density for the coefficient of the magnitude of the law effect for each state separately.

Also we will plot the actual log(homicide rate) together with the state specific counterfactual estimates.

lawco = group(s,"statelawcoef")

statenames = Dict()
for (st,stusab) in Iterators.zip(lawdates.STATE,lawdates.STUSAB)
    statenames[st] = stusab
    end

function plotlawts(s,lawdates,modeljoined)
let pl = []
    lawco = group(s,"statelawcoef")

    for i in @select(@subset(lawdates,:year .< 2022),:STATE).STATE
        den = density(s[:,Symbol("meanlawcoef"),:] .+ lawco[:,Symbol("statelawcoef[$i]"),:], title = "State $i = $(statenames[i])",
            legend=false,xlim=(-log(2),log(2)))
        plot!([(0.0,0.0),(0.0,5.0)], color="red",ylim=(0.0,5.0))
        push!(pl,den)
    end
    println("State law coefficient values")
    display(plot(pl...; size = (1000,1000)))

    display(density(s[:,:lawrate,:], title= "Law effect onset rate (1/yr)"))

    modeljoined.logdrate = log.(modeljoined.crudedeathrate)

    regions = Dict()
    for r in eachrow(stategroupdf)
        regions[r.STATE] = (group=r.stgroup,code=r.STUSAB)
    end
    pl = []
    samps = sample(1:500,10)

    for st in unique(modeljoined.STATE)
        sub = @subset(modeljoined, :STATE .== st)
        region = regions[st].group
        p = plot(sub.year,sub.logdrate,linewidth=3,title="$(st) = $(regions[st].code)",ylim=(-0.5,2.75))
        push!(pl,p)
        for samp in samps
            statelawcoef = s[samp,Symbol("statelawcoef[$st]"),1]
            statecoefs = [s[samp,Symbol("statecoefs[$st][$i]"),1] for i in 1:length(centers)]
            statebase = s[samp,Symbol("statebase[$st]"),1]
            regioncoefs = [s[samp,Symbol("regioncoefs[$region][$i]"),1] for i in 1:length(centers)]        
            
            lawrate = s[samp,:lawrate,1]
            startlaw = lawdates[lawdates.STATE .== st,:].year[1]
            startlaw = ismissing(startlaw) ? 3000 : startlaw
            years = 1979:2020
            pred = [statebase + 
                timeser(yr,regioncoefs .+ (statecoefpert .* statecoefs),centers,width) for yr in years]
            plot!(years,pred; color="orange",alpha=0.5)
            if startlaw < 3000
                plot!([(startlaw,-10.0),(startlaw,10.0)],color="red")
            end
            regpred = [statebase + 
                timeser(yr,regioncoefs ,centers,width) for yr in years]

        end
    end
    println("Each state shown with its locally estimated counterfactual (orange), actual (dark blue), and date of CC law if any (red vertical)")
    display(plot(pl...; size =(2000,3000), legend=false))
end

density(s[:,:meanlawcoef,:],title="Mean law coefficient")
end

plotlawts(s,lawdates,modeljoined)

State law coefficient values
Each state shown with its locally estimated counterfactual (orange), actual
 (dark blue), and date of CC law if any (red vertical)

Shall-issue vs permitless

One plausible explanation for the relative lack of clear measurable effect from permitless carry laws is that permitless actually may have had very little effect on how many people carry firearms in public. Many of these states were shall-issue and many people who wanted to carry may have had CCW permits issued before the permitless law passed. Therefore the change in number of people carrying may have been minimal.

Using the date on which states converted to shall-issue permits, with a similar analysis, we can determine whether more of an effect was visible after that legal change. Some states converted from no-issue to shall issue, while others converted from may-issue to shall-issue.

# dates collected from https://en.wikipedia.org/wiki/History_of_concealed_carry_in_the_United_States
shallissdates = CSV.read("data/ShallIssueDates.csv",DataFrame)

shallissue = leftjoin(fipscodes,shallissdates; on = :STUSAB)

silawdates = @chain leftjoin(DataFrame(STATE=1:55),shallissue; on=:STATE) begin
    @subset(:STATE .< 60)
    @rtransform(:year = ismissing(:ShallIssueDate) ? 3000 : :ShallIssueDate)
    @orderby(:STATE)
end

modl2 = guns(modeljoined.StateCode,maximum(modeljoined.StateCode),modeljoined.stgroup,maximum(modeljoined.stgroup),
        centers, width, modeljoined.year,log.(modeljoined.crudedeathrate),silawdates.year,statecoefpert)



savedfile = "./saved/shallisssamples.dat"
global sisam = []
if stat(savedfile).mtime > stat("./model.jl").mtime
    global sisam = deserialize(savedfile)
else
    global sisam = sample(modl2,NUTS(500,0.8),MCMCThreads(),500,3)
    serialize(savedfile,sisam)
end

Shall Issue coefficients

plotlawts(sisam, silawdates,modeljoined)

State law coefficient values
Each state shown with its locally estimated counterfactual (orange), actual
 (dark blue), and date of CC law if any (red vertical)

The Take Home

Permitless/Constitutional Carry laws

Many states that have passed constitutional carry laws have done so recently enough that there is no data on homicides available from CDC Wonder yet.

For those which passed the law long enough ago to have some post-law data, still many of them were in the last few years leaving only a few years of post-law data. The last few years is a time period where homicides have increased nearly nationwide and it is challenging to estimate the effect of a law when it is imposed on top of a nonlinearly changing trend. Estimates of a law effect created by comparing the actual data to counterfactuals estimated as small perturbations to the regional trend lead to estimates of the average law coefficient that has uncertain sign, with a most likely value of about 0.05 or so, meaning that gun homicides may increase by a factor of around 1.05, however the ability to resolve this number is so poor that the number could be anything from about -0.1 to +0.25 (multiplying the crime rate by anywhere from 0.9 to 1.28.

Any kind of convincing evidence for a consistent average net positive or negative effect of passing constitutional carry laws simply isn't there. What's more, we have decent bounds on the size of this effect for most states, it should be somewhere between perhaps -0.25 and 0.25 on the log scale, meaning approximately the rate of homicides would be multiplied by some number in the range 0.78 to 1.28.

Shall issue laws

Because many shall-issue laws were passed in the 1990's there is much better evidence for their long term effect on crime. However, using the same methodology as above the average effect is essentially symmetrically distributed around 0.0 with a posterior credible range of about -0.1 to 0.1. The only reasonably strong evidence we have in this data for a real effect is that Texas may have reduced crime with a coefficient around -0.3 and perhaps Pennsylvania caused an increase around 0.2, with Wisconsin possibly around 0.2. The PA estimate does see the counterfactual diverge around the time of the law, but the Wisconsin counterfactual diverges 20 years before the law, suggesting that that effect is really an artifact of a poor counterfactual estimate.

Firearms Sales

Firearms sales are driven by multiple factors. One is the fear of "bans" which drove a lot of sales around the time of the Sandy Hook school shooting in Dec 2012. And in general sales have increased over the last decade including a sudden increase after the turmoil of 2020, with women and racial minorities making up a historically unprecedented fraction of new sales (close to 50%).

Firearms sales can plausibly be both a cause of firearms crime (when people purchase firearms to commit crime), a response to violence (when people purchase firearms for defense after periods of crime), and even in some cases a deterrent to crime (when criminals meet with greater defensive responses than before). It is hard to say which of these effects dominates, and in fact most likely all of these occur at different times and different places. One thing is clear however, and that is that it is entirely possible for gun ownership to increase while crime decreases, or for crime to increase while gun ownership per capita decreases. Those who somehow wish to draw a direct causal link in the US between more guns and more gun crime simply do not have an easy obvious convincing argument.

The evidence suggests that in the post NYSRPA v Bruen era, forcing CA and NY and others to shall-issue their permits will by itself have very little effect on overall gun crime. More likely if there is an increase it will be due to the turmoil undergoing the whole country at the moment, with high levels of inflation and the war in Ukraine driving food prices up and crime in general increasing over the last few years.

Published from analyze.jmd using Weave.jl v0.10.10 on 2022-07-16.